### Replication code for 
### Eggers, A.E., Vivyan, N & Wagner, M. 'Corruption, accountability and gender: do female politicians face higher standards in public life?'
### The Journal of Politics

### (Written using R 3.2.4)


### set working directory and options

rm(list = ls()); gc()
setwd("~/Dropbox/Gender and corruption/replication_materials")


### load libraries

library(ggplot2)
library(rms)
library(stargazer)
library(tables)
library(lmtest)
library(multiwayvcov)
library(mvtnorm)


### Useful functions

mean2pct <- function(x){100*mean(x)}

meanlohi <- function(x){c(mean(x), quantile(x, probs = c(0.025, 0.975)))}

predict.rob <- function(object, vcov,newdata){
  tt <- terms(object)
  if(missing(newdata)){ newdata <- x$model }
  else {
    Terms <- delete.response(tt)
    m <- model.frame(Terms, newdata, na.action = na.pass, 
                     xlev = object$xlevels)
    if (!is.null(cl <- attr(Terms, "dataClasses"))) 
      .checkMFClasses(cl, m)
    X <- model.matrix(Terms, m, contrasts.arg = object$contrasts)
    offset <- rep(0, nrow(X))
    if (!is.null(off.num <- attr(tt, "offset"))) 
      for (i in off.num) offset <- offset + eval(attr(tt, 
                                                      "variables")[[i + 1]], newdata)
    if (!is.null(object$call$offset)) 
      offset <- offset + eval(object$call$offset, newdata)
    mmDone <- FALSE
  }
  m.coef <- object$coef
  fit <- as.vector(X %*% object$coef)
  se.fit <- sqrt(diag(X%*%vcov%*%t(X)))
  return(list(fit=fit,se.fit=se.fit))
} # Adapted from: https://stackoverflow.com/questions/3790116/using-clustered-covariance-matrix-in-predict-lm




########################


### Data preparation

## Load data
long.dat <- readRDS("experiment_data.rds")

## Create subset of data with observations from respondents' first choice task only
firstdat <- subset(long.dat, comparison == 1)



#########################


### Descriptive results for Appendix B

tab <- tabular( (Factor(resp.female, name = "Resp. gender",  levelnames = c("Male", "Female")))*(Factor(mp.female, name = "MP gender",  levelnames = c("Male", "Female")))*(Factor(mp.misconduct, name = "MP conduct",  levelnames = c("Good", "Bad"))) ~ 
                  Format(digits=3)*(Incumbent = voteinc)*(Share=mean2pct) +
                  (N=1) , data=long.dat)

for(outtype in c("html", "latex")){
  filename <- "tableB1"
  if(outtype == "html"){
    html(tab, paste0(filename, ".html"))
  } else {
    latex(tab, paste0(filename, ".tex"))
  }
}



#########################


### Regression analysis for Tables 1, C.1, E.1, and E.2

m1 <- lm(voteinc ~ mp.misconduct + mp.female + resp.female, long.dat)
se1 <- sqrt(diag(cluster.vcov(m1, long.dat$id))) # Cluster SEs

m1first <- update(m1, data = firstdat) # re-run on data from first choice-task only
se1first <- sqrt(diag(vcov(m1first, "HC3")))


m2 <- lm(voteinc ~ mp.misconduct + mp.female  + resp.female + 
           incid + chaid + minorid + socialgrade + agegroup +
           seat.type + factor(mp.age) + mp.prevjob + 
           factor(cha.age) + cha.prevjob + cha.sex
         , long.dat)
se2 <- sqrt(diag(cluster.vcov(m2, long.dat$id)))

m2first <- update(m2, data = firstdat) # re-run on data from first choice-task only
se2first <- sqrt(diag(vcov(m2first, "HC3")))


m3 <- lm(voteinc ~ mp.misconduct*resp.female + mp.female, long.dat)
se3 <- sqrt(diag(cluster.vcov(m3, long.dat$id)))

m3first <- update(m3, data = firstdat) # re-run on data from first choice-task only
se3first <- sqrt(diag(vcov(m3first, "HC3")))


m4 <- lm(voteinc ~ mp.misconduct*resp.female + mp.female + incid + chaid + minorid + socialgrade + agegroup+
           seat.type + factor(mp.age) + mp.prevjob + 
           factor(cha.age) + cha.prevjob + cha.sex
         , long.dat)
se4 <- sqrt(diag(cluster.vcov(m4, long.dat$id)))

m4first <- update(m4, data = firstdat) # re-run on data from first choice-task only
se4first <- sqrt(diag(vcov(m4first, "HC3")))


m5 <- lm(voteinc ~ mp.misconduct*mp.female + resp.female, long.dat)
se5 <- sqrt(diag(cluster.vcov(m5, long.dat$id)))

m5first <- update(m5, data = firstdat) # re-run on data from first choice-task only
se5first <- sqrt(diag(vcov(m5first, "HC3")))


m6 <- lm(voteinc ~ mp.misconduct*mp.female  + resp.female + incid + chaid + minorid + socialgrade + agegroup+
           seat.type + factor(mp.age) + mp.prevjob + 
           factor(cha.age) + cha.prevjob + cha.sex
         , long.dat)
se6 <- sqrt(diag(cluster.vcov(m6, long.dat$id)))

m6first <- update(m6, data = firstdat) # re-run on data from first choice-task only
se6first <- sqrt(diag(vcov(m6first, "HC3")))


m7 <- lm(voteinc ~ mp.misconduct*resp.female*mp.female, long.dat)
se7 <- sqrt(diag(cluster.vcov(m7, long.dat$id)))

m7first <- update(m7, data = firstdat) # re-run on data from first choice-task only
se7first <- sqrt(diag(vcov(m7first, "HC3")))


m8 <- lm(voteinc ~ mp.misconduct*resp.female*mp.female + incid + chaid + minorid + socialgrade + agegroup +
           seat.type + factor(factor(mp.age)) + mp.prevjob + 
           factor(cha.age) + cha.prevjob + cha.sex
         , long.dat)
se8 <- sqrt(diag(cluster.vcov(m8, long.dat$id)))

m8first <- update(m8, data = firstdat) # re-run on data from first choice-task only
se8first <- sqrt(diag(vcov(m8first, "HC3")))


## Make Table 1: Probability of voting for incumbent, MP misconduct and gender

for(outtype in c("text", "latex")){
  the.extension <- ifelse(outtype == "text", ".txt", ".tex")
  filename <- paste0("table1", the.extension) 
  stargazer(mget(paste0("m",c(1,5,3,7))),
            se = mget(paste0("se",c(1,5,3,7))),
            type = outtype,
            float = FALSE,
            out = filename,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            title = "", 
            omit.table.layout = "n",
            keep.stat = c("n", "rsq", "adj.rsq"),
            omit = "incid|chaid|minorid|socialgrade|agegroup|seat.type|mp.age|mp.prevjob|cha.age|cha.prevjob|cha.sex"
            , order = c("Constant", "^mp.misconduct$", "^mp.female$", "^resp.female$", 
                        "^mp.misconduct:mp.female$",
                        "^mp.misconduct:resp.female$",
                        "^resp.female:mp.female$"
            ),
            covariate.labels = c("Intercept", "MP misconduct", "MP female",
                                 "Respondent female", 
                                 "MP misconduct $\\times$ MP female",
                                 "MP misconduct $\\times$ Resp. female",
                                 "MP female $\\times$ Resp. female",
                                 "MP misconduct $\\times$ MP female $\\times$ Resp. female")
  )
  
}


## Make Table C.1: Probability of voting for incumbent MP by misconduct and gender, including models without controls

for(outtype in c("text", "latex")){
  the.extension <- ifelse(outtype == "text", ".txt", ".tex")
  filename <- paste0("tableC1", the.extension)
  
  stargazer(mget(paste0("m",c(1,2,5,6,3,4,7,8))),
            se = mget(paste0("se",c(1,2,5,6,3,4,7,8))),
            type = outtype,
            float = FALSE,
            out = filename,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            title = "", 
            omit.table.layout = "n",
            keep.stat = c("n", "rsq", "adj.rsq"),
            omit = "incid|chaid|minorid|socialgrade|agegroup|seat.type|mp.age|mp.prevjob|cha.age|cha.prevjob|cha.sex"
            , order = c("Constant", "^mp.misconduct$", "^mp.female$", "^resp.female$", 
                        "^mp.misconduct:mp.female$",
                        "^mp.misconduct:resp.female$",
                        "^resp.female:mp.female$"
            ),
            add.lines = list(c("Controls for voter, MP and chall. characteristics?", 
                               rep(c("No","Yes"), 4))),
            covariate.labels = c("Intercept", "MP misconduct", "MP female",
                                 "Respondent female", 
                                 "MP misconduct $\\times$ MP female",
                                 "MP misconduct $\\times$ Resp. female",
                                 "MP female $\\times$ Resp. female",
                                 "MP misconduct $\\times$ MP female $\\times$ Resp. female")
  )
  
}


## Table E.1: Robustness of main regression results in first choice task subsample

for(outtype in c("text", "latex")){
  the.extension <- ifelse(outtype == "text", ".txt", ".tex")
  filename <- paste0("tableE1", the.extension)
  
  stargazer(mget(paste0(paste0("m",rep(c(1,5,3,7),each = 2)),c("", "first"))),
            se = mget(paste0(paste0("se",rep(c(1,5,3,7),each = 2)),c("", "first"))),
            type = outtype,
            float = FALSE,
            out = filename,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            column.labels = rep(c("All tasks", "First task"), 4),
            title = "", 
            omit.table.layout = "n",
            keep.stat = c("n", "rsq", "adj.rsq"),
            omit = "incid|chaid|minorid|socialgrade|agegroup|seat.type|mp.age|mp.prevjob|cha.age|cha.prevjob|cha.sex"
            , order = c("Constant", "^mp.misconduct$", "^mp.female$", "^resp.female$", 
                        "^mp.misconduct:mp.female$",
                        "^mp.misconduct:resp.female$",
                        "^resp.female:mp.female$"
            ),
            covariate.labels = c("Intercept", "MP misconduct", "MP female",
                                 "Respondent female", 
                                 "MP misconduct $\\times$ MP female",
                                 "MP misconduct $\\times$ Resp. female",
                                 "MP female $\\times$ Resp. female",
                                 "MP misconduct $\\times$ MP female $\\times$ Resp. female")
  )
  
}



## Table E.2: Robustness of main regression results in first choice task subsample, incorporating controls 

for(outtype in c("text", "latex")){
  the.extension <- ifelse(outtype == "text", ".txt", ".tex")
  filename <- paste0("tableE2", the.extension)
  
  stargazer(mget(paste0(paste0("m",rep(c(2,6,4,8),each = 2)),c("", "first"))),
            se = mget(paste0(paste0("se",rep(c(2,6,4,8),each = 2)),c("", "first"))),
            type = outtype,
            float = FALSE,
            out = filename,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            column.labels = rep(c("All tasks", "First task"), 4),
            title = "", 
            omit.table.layout = "n",
            keep.stat = c("n", "rsq", "adj.rsq"),
            omit = "incid|chaid|minorid|socialgrade|agegroup|seat.type|mp.age|mp.prevjob|cha.age|cha.prevjob|cha.sex"
            , order = c("Constant", "^mp.misconduct$", "^mp.female$", "^resp.female$", 
                        "^mp.misconduct:mp.female$",
                        "^mp.misconduct:resp.female$",
                        "^resp.female:mp.female$"
            ),
            add.lines = list(c("Controls for voter, MP and chall. characteristics?", 
                               rep(c("Yes","Yes"), 4))),
            covariate.labels = c("Intercept", "MP misconduct", "MP female",
                                 "Respondent female", 
                                 "MP misconduct $\\times$ MP female",
                                 "MP misconduct $\\times$ Resp. female",
                                 "MP female $\\times$ Resp. female",
                                 "MP misconduct $\\times$ MP female $\\times$ Resp. female")
  )
  
}




### Regression analysis for Table D.1

m4extra <- lm(voteinc ~ mp.misconduct*resp.female + mp.female + 
                mp.misconduct*incid + 
                mp.misconduct*chaid + 
                mp.misconduct*minorid + 
                mp.misconduct*socialgrade + 
                mp.misconduct*agegroup+
           seat.type + factor(mp.age) + mp.prevjob + 
           factor(cha.age) + cha.prevjob + cha.sex
         , long.dat)
se4extra <- sqrt(diag(cluster.vcov(m4extra, long.dat$id)))
coeftest(m4extra, cluster.vcov(m4extra, long.dat$id))


m8extra <- lm(voteinc ~ mp.misconduct*resp.female*mp.female + 
                mp.misconduct*resp.female*incid + 
                mp.misconduct*resp.female*chaid + 
                mp.misconduct*resp.female*minorid + 
                mp.misconduct*resp.female*socialgrade + 
                mp.misconduct*resp.female*agegroup +
                seat.type + factor(factor(mp.age)) + mp.prevjob + 
                factor(cha.age) + cha.prevjob + cha.sex
              , long.dat) # model with interactions with all controls
se8extra <- sqrt(diag(cluster.vcov(m8extra, long.dat$id)))
coeftest(m8extra, cluster.vcov(m8extra, long.dat$id))


## Make Table D.1: Probability of voting for incumbent MP by misconduct and gender, controlling for voter attribute $\times$ treatment interactions

for(outtype in c("text", "latex")){
  the.extension <- ifelse(outtype == "text", ".txt", ".tex")
  filename <- paste0("tableD1", the.extension)
  
  stargazer(mget(paste0(paste0("m",rep(c(4,8),each = 2)),c("", "extra"))),
            se = mget(paste0(paste0("se",rep(c(4,8),each = 2)),c("", "extra"))),
            type = outtype,
            float = FALSE,
            out = filename,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            title = "", 
            omit.table.layout = "n",
            keep.stat = c("n", "rsq", "adj.rsq"),
            omit = "incid|chaid|minorid|socialgrade|agegroup|seat.type|mp.age|mp.prevjob|cha.age|cha.prevjob|cha.sex|^mp.misconduct$|^mp.female$|^resp.female$|^mp.misconduct:mp.female$|^resp.female:mp.female$|Constant"
            , order = c(
              "^mp.misconduct:resp.female$"
            ),
            add.lines = list(c("Controls for main effects of voter, MP and chall. characteristics?", 
                               rep(c("Yes","Yes"), 4)), 
                             c("Controls for voter characteristic $\\times$ MP misconduct interactions?", 
                               c("No", "Yes", "No", "Yes")),
                             c("Controls for voter characteristic $\\times$ MP misconduct $\\times$ MP gender interactions and constituent terms?", 
                               c("No", "No", "No", "Yes"))
            )
            ,covariate.labels = c("MP misconduct $\\times$ Resp. female",
                                  "MP misconduct $\\times$ MP female $\\times$ Resp. female")
  )
}



### Regression analysis for Table F.1

the.dat <- subset(long.dat, resp.female == 0)
m10 <- lm(voteinc ~ mp.misconduct*mp.female*cha.female + incid + chaid + minorid + socialgrade + agegroup +
           seat.type + factor(mp.age) + mp.prevjob + 
           factor(cha.age) + cha.prevjob
         , the.dat)
se10 <- sqrt(diag(cluster.vcov(m10, the.dat$id)))

m10nocontrol <- lm(voteinc ~ mp.misconduct*mp.female*cha.female 
          , the.dat)
se10nocontrol <- sqrt(diag(cluster.vcov(m10nocontrol, the.dat$id)))

the.dat <- subset(long.dat, resp.female == 1)
m11 <- lm(voteinc ~ mp.misconduct*mp.female*cha.female + incid + chaid + minorid + socialgrade + agegroup+
           seat.type + factor(mp.age) + mp.prevjob + 
           factor(cha.age) + cha.prevjob
         , the.dat)
se11 <- sqrt(diag(cluster.vcov(m11, the.dat$id)))

m11nocontrol <- lm(voteinc ~ mp.misconduct*mp.female*cha.female 
          , the.dat)
se11nocontrol <- sqrt(diag(cluster.vcov(m11nocontrol, the.dat$id)))


## Make Table F1: Probability of voting for incumbent, MP misconduct, MP gender and challenger gender

for(outtype in c("text", "latex")){
  the.extension <- ifelse(outtype == "text", ".txt", ".tex")
  filename <- paste0("tableF1", the.extension)
  stargazer(mget(paste0("m",c("10nocontrol","11nocontrol", "10", "11"))),
            se = mget(paste0("se",c("10nocontrol","11nocontrol", "10",  "11"))),
            type = outtype,
            float = FALSE,
            out = filename,
            dep.var.caption = "",
            dep.var.labels.include = FALSE,
            column.labels = rep(c("Male resp.", "Female resp."),2),
            column.separate = c(1,1,1,1),
            title = "", 
            omit.table.layout = "n",
            intercept.bottom = FALSE, intercept.top = TRUE,
            keep.stat = c("n", "rsq", "adj.rsq"),
            omit = "incid|chaid|minorid|socialgrade|agegroup|seat.type|mp.age|mp.prevjob|cha.age|cha.prevjob",
            add.lines = list(c("Controls for voter, MP and chall. characteristics?", 
                               rep(c("No", "Yes"), each = 2))),
            order = c("Constant", "^mp.misconduct$", "^mp.female$", 
                      "^cha.female$","^resp.female$",
                      paste0("^mp.misconduct:", c("mp.female$","cha.female$", "resp.female$"))),
                      covariate.labels = c("Intercept", "MP misconduct", "MP female",
                                 "Challenger female", 
                                 "MP misconduct $\\times$ MP female",
                                 "MP misconduct $\\times$ Chall. female",
                                 "MP female $\\times$ Chall. female",
                                 "MP misconduct $\\times$ MP female $\\times$ Chall. female")
  )
}





#########################

### Make Figures 2 and F.1, based above regression models


## Make Figure 2a showing treatment effects of misconduct by MP gender and voter gender

# do this by simulation

m8.betasim <- rmvnorm(n = 10000, mean = coef(m7), sigma = cluster.vcov(m7, long.dat$id))

# now get marginal effect of mp.misconduct by respondent and MP gender
plotdf <- data.frame(expand.grid(mp.female = c(0,1,"diff"), resp.female = c(0,1)),
                     eff = NA, lo = NA, hi = NA)

plotdf[plotdf$mp.female == 0 & plotdf$resp.female == 0, 3:5] <- meanlohi(m8.betasim[,"mp.misconduct"])
plotdf[plotdf$mp.female == 1 & plotdf$resp.female == 0, 3:5] <- meanlohi(
  rowSums(m8.betasim[,c("mp.misconduct", "mp.misconduct:mp.female")]))
plotdf[plotdf$mp.female == "diff" & plotdf$resp.female == 0, 3:5] <- meanlohi(
  rowSums(m8.betasim[,c("mp.misconduct", "mp.misconduct:mp.female")]) - m8.betasim[,"mp.misconduct"]
  )

plotdf[plotdf$mp.female == 0 & plotdf$resp.female == 1, 3:5] <- meanlohi(
  rowSums(m8.betasim[,c("mp.misconduct", "mp.misconduct:resp.female")]))
plotdf[plotdf$mp.female == 1 & plotdf$resp.female == 1, 3:5] <- meanlohi(
  rowSums(m8.betasim[,c("mp.misconduct", "mp.misconduct:mp.female", "mp.misconduct:resp.female", "mp.misconduct:resp.female:mp.female")]))
plotdf[plotdf$mp.female == "diff" & plotdf$resp.female == 1, 3:5] <- meanlohi(
  rowSums(m8.betasim[,c("mp.misconduct", "mp.misconduct:mp.female", "mp.misconduct:resp.female", "mp.misconduct:resp.female:mp.female")]) - rowSums(m8.betasim[,c("mp.misconduct", "mp.misconduct:resp.female")])
  )


plotdf$mp.female.neat <- car:::recode(plotdf$mp.female, '"0" = "Male MP"; "1" = "Female MP"; 
                                      "diff" = "Difference in effects\\n(Female MP - Male MP)"',
                                   as.factor.result = TRUE, 
                                   levels = c("Difference in effects\n(Female MP - Male MP)", "Female MP", "Male MP"))
plotdf$resp.female.neat <- car:::recode(plotdf$resp.female, '0 = "Male voter"; 1 = "Female voter"',
                                   as.factor.result = TRUE, 
                                   levels = c("Male voter", "Female voter"))

ggplot(plotdf, aes(x = mp.female.neat, y = eff)) + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.5) +
  geom_linerange(aes(x=mp.female.neat, xend=mp.female.neat, ymin=lo, ymax=hi), size = 0.6) +
  geom_point(size = 3.5, shape = 21, fill = "white") +
  labs(x = "", y = "Effect of MP misconduct on support for incumbent") + 
  coord_flip() +
  facet_wrap( ~ resp.female.neat, ncol = 1) + 
  theme_bw()
ggsave("figure2a.pdf", width = 5, height = 3.5)



## Make Figure 2b showing predicted support by MP conduct, MP gender and voter gender

newdf <- data.frame(expand.grid(mp.misconduct = c(0,1), resp.female = c(0,1),
                                mp.female = c(0,1)),
                    incid = 1, chaid = 0, minorid = 0, 
                    socialgrade = factor("ABC1", levels = levels(long.dat$socialgrade)),
                    agegroup = factor("40-59", levels = levels(long.dat$agegroup)),
                    seat.type = factor("ConLab",levels = levels(long.dat$seat.type)),
                    mp.age = factor("45", levels = levels(factor(long.dat$mp.age))),
                    mp.prevjob = factor("a business manager", levels = levels(long.dat$mp.prevjob)),
                    cha.age = factor("51", levels = levels(factor(long.dat$cha.age))),
                    cha.prevjob = factor("a teacher", levels = levels(long.dat$cha.prevjob)),
                    cha.sex = factor("Male", levels = levels(long.dat$cha.sex)),
                    voteinc = c(0,1)
)
preds <- predict.rob(m7, vcov = cluster.vcov(m7, long.dat$id), newdata = newdf)


newdf$yhat <- preds$fit
newdf$se.yhat <- preds$se.fit
newdf$lo <- newdf$yhat - (1.96*newdf$se.yhat)
newdf$hi <- newdf$yhat + (1.96*newdf$se.yhat)
newdf$mp.misconduct.neat <- car:::recode(newdf$mp.misconduct, '0 = "Good conduct"; 1 = "Misconduct"',
                                         as.factor.result = TRUE, 
                                         levels = c("Good conduct", "Misconduct"))
newdf$mp.female.neat <- car:::recode(newdf$mp.female, '0 = "Male MP"; 1 = "Female MP"',
                                     as.factor.result = TRUE, 
                                     levels = c("Female MP","Male MP"))
newdf$resp.female.neat <- car:::recode(newdf$resp.female, '0 = "Male voter"; 1 = "Female voter"',
                                       as.factor.result = TRUE, 
                                       levels = c("Male voter", "Female voter"))

ggplot(newdf, aes(x = mp.female.neat, y = yhat, fill = mp.misconduct.neat)) + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.5) +
  geom_linerange(aes(x=mp.female.neat, xend=mp.female.neat, ymin=lo, ymax=hi), size = 0.6,
                 position = position_dodge(width = - 0.7), show_guide=F) +
  geom_point(position = position_dodge(width = - 0.7), size = 3.5, shape = 21, show_guide = TRUE) +
  labs(x = "", y = "Predicted support for incumbent") + 
  coord_flip() +
  facet_wrap( ~ resp.female.neat, ncol = 1) + 
  theme_bw() +
  ylim(0.1,0.5) +
  scale_fill_manual(name = "MP conduct",  
                    values = c("White","Black")) +
  theme(legend.position = "bottom")
ggsave("figure2b.pdf", width = 5, height = 3.5)



## Make Figure F.1 showing treatment effects of misconduct by voter, incumbent and challenger gender

# do this by simulation

the.maledat <- subset(long.dat, resp.female == 0)
the.femaledat <- subset(long.dat, resp.female == 1)

m10.betasim <- rmvnorm(n = 10000, mean = coef(m10), sigma = cluster.vcov(m10, the.maledat$id)) # model for male resp
m11.betasim <- rmvnorm(n = 10000, mean = coef(m11), sigma = cluster.vcov(m11, the.femaledat$id)) # model for female resp

# now get marginal effect of mp.misconduct by respondent and MP gender
plotdf <- data.frame(expand.grid(mp.female = c(0,1,"diff"), cha.female = c(0,1), resp.female = c(0,1)),
                     eff = NA, lo = NA, hi = NA)

meanlohi <- function(x){c(mean(x), quantile(x, probs = c(0.025, 0.975)))}

# male respondents
plotdf[plotdf$mp.female == 0 & plotdf$cha.female == 0 & plotdf$resp.female == 0, 4:6] <- meanlohi(m10.betasim[,"mp.misconduct"])
plotdf[plotdf$mp.female == 1 & plotdf$cha.female == 0 & plotdf$resp.female == 0, 4:6] <- meanlohi(
  rowSums(m10.betasim[,c("mp.misconduct", "mp.misconduct:mp.female")])
  )
plotdf[plotdf$mp.female == "diff" & plotdf$cha.female == 0 & plotdf$resp.female == 0, 4:6] <- meanlohi(
  rowSums(m10.betasim[,c("mp.misconduct", "mp.misconduct:mp.female")]) - m10.betasim[,"mp.misconduct"]
)
plotdf[plotdf$mp.female == 0 & plotdf$cha.female == 1 & plotdf$resp.female == 0, 4:6] <- meanlohi(
  rowSums(m10.betasim[,c("mp.misconduct", "mp.misconduct:cha.female")])
)
plotdf[plotdf$mp.female == 1 & plotdf$cha.female == 1 & plotdf$resp.female == 0, 4:6] <- meanlohi(
  rowSums(m10.betasim[,c("mp.misconduct", "mp.misconduct:mp.female", "mp.misconduct:cha.female", "mp.misconduct:mp.female:cha.female")])
)
plotdf[plotdf$mp.female == "diff" & plotdf$cha.female == 1 & plotdf$resp.female == 0, 4:6] <- meanlohi(
  rowSums(m10.betasim[,c("mp.misconduct", "mp.misconduct:mp.female", "mp.misconduct:cha.female", "mp.misconduct:mp.female:cha.female")]) - rowSums(m10.betasim[,c("mp.misconduct", "mp.misconduct:cha.female")])
)

# female respondents
plotdf[plotdf$mp.female == 0 & plotdf$cha.female == 0 & plotdf$resp.female == 1, 4:6] <- meanlohi(m11.betasim[,"mp.misconduct"])
plotdf[plotdf$mp.female == 1 & plotdf$cha.female == 0 & plotdf$resp.female == 1, 4:6] <- meanlohi(
  rowSums(m11.betasim[,c("mp.misconduct", "mp.misconduct:mp.female")])
)
plotdf[plotdf$mp.female == "diff" & plotdf$cha.female == 0 & plotdf$resp.female == 1, 4:6] <- meanlohi(
  rowSums(m11.betasim[,c("mp.misconduct", "mp.misconduct:mp.female")]) - m11.betasim[,"mp.misconduct"]
)
plotdf[plotdf$mp.female == 0 & plotdf$cha.female == 1 & plotdf$resp.female == 1, 4:6] <- meanlohi(
  rowSums(m11.betasim[,c("mp.misconduct", "mp.misconduct:cha.female")])
)
plotdf[plotdf$mp.female == 1 & plotdf$cha.female == 1 & plotdf$resp.female == 1, 4:6] <- meanlohi(
  rowSums(m11.betasim[,c("mp.misconduct", "mp.misconduct:mp.female", "mp.misconduct:cha.female", "mp.misconduct:mp.female:cha.female")])
)
plotdf[plotdf$mp.female == "diff" & plotdf$cha.female == 1 & plotdf$resp.female == 1, 4:6] <- meanlohi(
  rowSums(m11.betasim[,c("mp.misconduct", "mp.misconduct:mp.female", "mp.misconduct:cha.female", "mp.misconduct:mp.female:cha.female")]) - rowSums(m11.betasim[,c("mp.misconduct", "mp.misconduct:cha.female")])
)


plotdf$mp.female.neat <- car:::recode(plotdf$mp.female, '"0" = "Male MP"; "1" = "Female MP"; 
                                      "diff" = "Difference in effects\\n(Female MP - Male MP)"',
                                      as.factor.result = TRUE, 
                                      levels = c("Difference in effects\n(Female MP - Male MP)", "Female MP", "Male MP"))
plotdf$cha.female.neat <- car:::recode(plotdf$cha.female, '0 = "Male challenger"; 1 = "Female challenger"',
                                       as.factor.result = TRUE, 
                                       levels = c("Male challenger", "Female challenger"))
plotdf$resp.female.neat <- car:::recode(plotdf$resp.female, '0 = "Male voter"; 1 = "Female voter"',
                                        as.factor.result = TRUE, 
                                        levels = c("Male voter", "Female voter"))

ggplot(plotdf, aes(x = mp.female.neat, y = eff, fill = cha.female.neat)) + 
  geom_hline(yintercept = 0, linetype = "dashed", size = 0.5) +
  geom_linerange(aes(x=mp.female.neat, xend=mp.female.neat, ymin=lo, ymax=hi), 
                 position = position_dodge(width = - 0.7),
                 show_guide = FALSE, size = 0.6) +
  geom_point(position = position_dodge(width = - 0.7), 
             show_guide = TRUE, size = 3.5, shape = 21) +
  labs(x = "", y = "Effect of MP misconduct on support for incumbent") + 
  coord_flip() +
  facet_wrap( ~ resp.female.neat, ncol = 1) + 
  theme_bw() +
  scale_fill_manual(name = "",  
                    values = c("Black","white")) +
  theme(legend.position = "bottom")
ggsave("figureF1.pdf", width = 7.5, height = 6)





